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The site-diluted transverse field Ising model in two dimensions is studied with Quantum- 
Monte-Carlo simulations. Its phase diagram is determined in the transverse field (T) and 
temperature (T) plane for various (fixed) concentrations (p). The nature of the quantum 
Griffiths phase at zero temperature is investigated by calculating the distribution of the local 
zero-frequency susceptibility. It is pointed out that the nature of the Griffiths phase is different 
for small and large I\ 
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§1. Introduction 

In numerous investigations of the effect of quantum fluctuations on the order-disorder phase 
transition the transverse-field spin-1/2 Ising model with ferromagnetic interactions has served as 
the simplest model on which these effects can be studied.0) This model is of particular interest, 
since the strength of the quantum fluctuations can be controlled by varying the strength T of the 
transverse field: a pure quantum phase transition (QPT) occurs at zero temperature when the 
external field exceeds a critical value. The necessary prerequisite for the existence of long range 
order is the existence of an infinite connected component or cluster in the lattice, because only an 
infinite system can have a spontaneous magnetization. As has been recognized by Harris in the 
context of randomly site diluted lattices an infinite connected cluster is not only necessary but also 
sufficient for the existence of such a QPT. Any percolating cluster is at least one-dimensional (i.e. 
the fractal dimension of its backbone is larger than one) . Since even the one-dimensional transverse 
field Ising model shows a QPT at some non-vanishing strength of the transverse field, any infinite 
cluster must have a QPT. 

The investigation of the effects of randomness on the transverse field Ising models (TFI) has 
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attracted a lot of interest recently]!' In particular, the peculiar features of random systems which 
derive from the Griffiths-McCoy singularity have been studied intensively.!^ 

As Griffiths has pointed out^ randomly diluted (classical) ferromagnets have a non-analytic 
free energy even away from the thermodynamical critical temperature for a given concentration p. 
As a matter of fact for any fixed concentration p of occupied sites such a non-analytic behavior 
persists up to a temperature which coincides with the transition temperature of the pure case (T c 
for p = 1). The reason for these anomalies is that there is a non- vanishing probability to find 
arbitrarily large clusters of connected spins for any concentration p. Since large clusters tend to 
order ferromagnetically below T c , they act very coherently. This causes singularities in the response 
to an external field and in the dynamical properties even above the critical point. These effects 
have been pointed out as significant effects in the distribution of zeros of the partition function and 
the Laplace transform of dynamical correlation .!) 

The Griffiths phase persists in the presence of a transverse field and is, for fixed concentration 
p, located between the ferromagnetic phase boundary T c (T;p) and the pure paramagnetic phase 
boundary T c (T;p = 1). While the effect of the Griffiths-McCoy singularity on the relaxation process 
of the diluted Ising model causes a non-exponential decay of the time correlation function, i^'&'i* 
it has no effect on the thermodynamic properties for non-zero temperatures. At zero temperature, 
however, the singularity causes singular behavior even in static propertieso'li-f 'li-iP' .Par- 
ticularly noteworthy are the divergence of the local susceptibility and the algebraic decay of the 
dynamical correlations. 

In this paper we study the transverse field spin-1/2 Ising model on a randomly diluted square 
lattice. Here, we mainly investigate the distribution of the local susceptibility using the continu- 
ous (imaginary) time cluster algorithm for transverse Ising models developed by Kawashima and 
Rieger0) (continuous time cluster algorithms for various other quantum mechanical models have 
been proposed recently^). It turns out that this distribution has a power law decay for large 
values, and consequently the local susceptibility diverges at some points in the disordered phase. 
A similar behavior has been found for quantum spin glass systems, but for these systems only the 
local nonlinear susceptibilities diverge .0)>lll|) The exponent for the power law decay of the distribu- 
tion is related to the dynamical exponent z. As our most important result we confirm the result of 
Senthil and SachdevEI) that at the quantum critical point time scales diverge exponentially fast at 
the percolation threshold, implying that the dynamical critical exponent z rr n is infinite, whereas in 
2 or 3-dinaensiooa, t_ L ising spin passes a nni.e vaje for J nas been estanatad.B 
The concentration dependence of the phase boundary in the site or bond diluted square lattices at 
zero temperature is also of interest. Particularly interesting features of the phase diagram such as 
the existence of a multi-critical point and a straight vertical phase boundary have been observed 
for small transverse fields. For such small transverse fields the percolation threshold de- 
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termines the boundary between the ordered (ferromagnetic) phase and the disordered phase. On 
the other hand for large transverse field, even in the geometrically connected cluster (i.e. above the 
percolation threshold), the magnetically connected cluster is reduced by quantum fluctuations. 

§2. Model and Method 

The model that we consider here is the spin- 1/2 Ising model in a transverse field on a square lattice 
with random site dilution, i.e. sites on a square lattice are occupied with spins with probability p 
and empty with probability 1 — p. Only occupied nearest neighbor sites interact. The model is 
thus defined by the quantum mechanical Hamiltonian 

h = -jJ2 - r £ etaf , (i) 

where a% are Pauli spin matrices, denotes nearest neighbor pairs on a L x L square lattice with 
periodic boundary conditions, J is the ferromagnetic interaction strength and T is the strength of 
the transverse field; £j G {0, 1} are quenched random variables modeling the dilution: £j = 1 with 
probability p and £j = with probability 1 — p. 

To map the 2-dimensional quantum system to the 3-dimensional classical system, we use the 
Suzuki- Trotter decomposition. Then the free energy T of the system (|l]) at inverse temperature 
(3 = 1/T is obtained as the limit of a 3-dimensional classical Ising model: 

T = -(3~ l lim In Tr exp(-S class ) , (2) 

n— >oo 

with 

n 

£iEjSj(T)Sj{T) - K ver y^£jgj(r)gj(r + 1) , (3) 

r=l (ij) T ,i 

and 

Khor = At J , 

K ver = -ilntanh(Arr) , (4) 
At = 0/n . 

The classical action (^|) is the Hamiltonian of a cubic lattice of L x L x n classical Ising spins 
Si{r) = ±1. The additional index r = l,...,n labels the n two-dimensional (imaginary) time 
slices within which spins interact via Kh or and among which they interact with strength K ver . The 
number of time slices n is called the Trotter number and the third (imaginary time) axis is called 
the Trotter axis. 

Note that an empty (occupied) site i in one two-dimensional time slice implies a whole column 
of empty (occupied) sites in the Trotter direction: the quenched disorder is perfectly correlated 
in the imaginary time direction. Therefore the model we study is not a diluted ferromagnet on 
a cubic lattice with uncorrelated disorder, and its universality class will be different from that of 
the site diluted cubic ferromagnet. We obtain the phase diagram for the diluted model by the 
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conventional quantum Monte Carlo method using the world line update for finite Trotter numbers, 
Fig. 1. The phase boundaries for p = 0.8, 0.7 and 0.6 are also shown. For finite temperatures, the 
phase boundary shrinks as p decreases towards p c (~ 0.59). It is shown that the Griffiths phase for 
fixed concentration p is located between the FM phase boundary T c (T;p) and the pure PM phase 
boundary T c (T;p = 1). But the conventional method has difficulties in the regions close to T = 0, 
because we need n proportional to [5 to take the Trotter limit explicitly. Moreover, it is extremely 
hard to equilibrate the system for weak transverse fields. In this work, we use the new continuous 
time cluster algorithm developed by Kawashima and Rieger0) to avoid these difficulties. For a 
description of the algorithm and details of the implementation we refer to.0) 

First, we have checked our method for the pure case (i.e. p = 1), for which the phase diagram 
is already known and the critical field value and the thermal exponents have been estimated from 
series expansions.©'© As a test of the cluster algorithm, we plot the square of the magnetization 
(M 2 ) obtained at each Monte Carlo step, Fig. 2. At high temperature and low T the system is 
in the disordered region. However, due to the strong coupling constant along the Trotter axis, the 
time to reach equilibration is large. The use of the cluster algorithm reduces this time significantly. 
In Fig. 3, we study n dependence of the critical values of V obtained by the Binder plot ( the 
dimensionless ratio of moments of the order parameter ) at [3 = 10. The value obtained from 
the continuous cluster algorithm is approximately 3.05. This is indeed close to the value which is 
obtained by extrapolating from the data of the conventional method. Finally, we follow Kawashima 
and Rieger@) and estimate the critical T at zero temperature and to obtain the critical exponents 
from finite size scaling. 

Close to the quantum critical point, quantities are expected to obey the finite size scaling form 



where 5 = T — T c and a is the finite size scaling exponent of the quantity A. From the equivalence 
with the 3-dimensional classical Ising model one knows the dynamical exponent to be z = 1. Thus, 
we can perform conventional one-parameter finite size scaling, if we choose the aspect ratio (3/L 
to be constant. We work with a constant aspect ratio (3/L = 1/5 and estimate the magnetization 
m (a = —f3/v), the uniform susceptibility x ( a = 2 — r/) and the Binder ratio g (a = 0) in order 
to determine the values of v, r/ and (5, see Fig. 4. As a result, we estimate r c ~ 3.06, v ~ 0.64, 
and f3 ~ 0.31 from m, r c ~ 3.07, v ~ 0.63 and r\ ~ 0.04 from \, and T c ~ 3.08 and v ~ 0.63 from 
g. These values agree well with the series expansion resultsil' and with those obtained earlier by 
Kawashima and RiegeJiil' using the same method. 

§3. Quantum Griffiths phase 

The phase boundary in the site diluted square lattices is depicted in Fig. 5. As has been 
mentioned above we expect the Griffiths phase in the classical model (r = 0) between the critical 




(5) 
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temperature of the pure model p = 1 and the thermodynamical critical temperature for the concen- 
tration T c (p). The singularity in this Griffiths phase is due to the finite probability of large clusters 
which behave coherently below T c . For low concentration the probability is given approximately 
by 

P(N)dN ocp N dN = exp(-A |lnp|) , (6) 

i.e. it vanishes exponentially with the number of spins. On the other hand, at T = (quantum 
region) the existence of a vertical line (parallel to the T-axis) at p c (the percolation threshold 
~ 0.590)) extending from r = to Tm follows from the nature of the backbone of the percolating 
cluster,© which has a fractal dimension between 1 and 2. Therefore the critical transverse field 
strength of the pure one-dimensional transverse Ising model is a lower bound for Tm, i-e. Tm > 
J = 1. 

Again, due to the presence of arbitrarily large connected clusters the entire area below T c is 
Griffiths phase (quantum Griffiths phase) even the ordered phase contains anomalies. The param- 
agnetic phase lies above T c . Since we are at zero temperature close to the quantum phase transition 
points T c (p), we encounter as a new feature much stronger singularities than in the classical case, 
a fact that has first been noted by McCoyS and has been intensely investigated recently in various 
situationsS©©'@)'§ The transition across the vertical line including the multi-critical point 
at (0, 1 — pc, Tm) as well as the Griffiths-McCoy singularity for small transverse fields T <Tm have 
also been discussed quite recently. 

The main aims of our present investigation are to determine the dynamical exponent in the 
quantum Griffiths phase and to confirm the conjecture of the straight vertical phase boundary. For 
P "C p c the probability distribution of the existence of a spherical dense cluster with TV spins is 
given by Eq.(|6|). In the transformed 3D classical Ising model, such a cluster forms a rod along the 
Trotter axis as shown in Fig. 6. The spins in the cluster are strongly coupled in real space and 
may have a domain wall perpendicular to the Trotter axis. The upward or downward solid arrow 
indicates parallel or anti-parallel with respect to the state of the cluster in real space. The insertion 
of a domain wall costs an energy AE 

AE ~ CN , (7) 

where C is stiffness constant. The probability for such a domain wall decreases exponentially as 
exp(— CN). Therefore the correlation length along the Trotter axis as a function of a cluster size 
A is given by 

£ T (iV) ~ exp (CN) . (8) 

Now we introduce the local susceptibility Xiocah which is the response of the expectation value of 
a spin af to a local (longitudinal) field hi on site i which is expressed by hiaf in the Hamiltonian. 

d 



X,local • (9) 
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In the continuous time method, following Kawashima and Rieger the expectation value for the 
local susceptibility is 

Xiocai= I K(rK(0))=/3(m?) . (10) 
Jo 

Here, the local magnetization rrtj at site i is the difference between the total length of all +segments 
and the total length of all —segments divided by (3. From the relation (|8|), xiocal is proportional to 
£ T . Now let us consider the distribution of £ T , P (xiocal)- The dependence of the correlation length 
^ T on the cluster size ./V is given by (|8|) and N has the distribution (||). Thus the distribution of £ T 
is 

P(t T )dt T ~&d£ T , (11) 

where 

* = 7?-i- <*> 

Consequently P (xiocal) nas a power law dependence on xiocal at large values, 

P (Xiocal) dXlocal - {Xlocalf dxiocal • (13) 

The numerical estimation of the integrated distribution function has a better statistics than 
P(xiocai)- Therefore, in the practical calculation we obtain the integrated distribution function 

F {Xiocal) = P(t)dt~±- (Xiocal)^ ' ■ (14) 

J Xiocal m P 

In Fig. 7-11, the integrated distribution functions are plotted for various values of T, and p. We 
expect that the integrated distribution of Xiocal will be cut off at Xiocal = P, and that increasing f3 
will simply extend the range over which the data fall on a straight line. 
Near the percolation threshold p c , where the form (||) is no longer valid, P (N) is given by© 

P(iV)~iV-ex P (-^) . (15) 
Here, £ is correlation length in real space which diverges as £ ~ \p — p c \~ l ' p , and D is the fractal 



dimension of the percolating cluster in two space dimensions. Using Eq.(15) instead of Eq.@, one 
obtains a formula for F (xiocal) that is valid near the percolation threshold as long as the phase 
boundary is vertical, i.e. T < Tm- 

F (Xiocal) ~— (Xiocal)-^ . (16) 

a 

Following Rieger and YoungS one can relate the power of the tail of F(xiocai) to a single dynamical 
exponent z(p, T) that varies continuously with p and T, 

In Fixiocal ) = - - In xiocal + const . , (17) 

z 



6 



where d is the space dimension, i.e. d=2 here. From fll4|), ( Jl6| ) and (O) it follows that 

d — \np/C for p ^ p c 



z 



(18) 

a/C£ for p — > p c . 



The second relation implies that 2 diverges at p c algebraically© like 

2 oc \p — p c \~ DUp for p — > p c and r < Tm- (19) 

In Fig. 11 the data for d/z are plotted for various values of T, and p. For T = 0.7, 1.0 (< Tm), 
the data for d/z vanish at p c , whereas for T = 1.5,2.0, they vanish for p > p c . For T > Tm the 
phase boundary is no longer parallel to the T-axis, and the quantum fluctuations are strong enough 
to destroy the ferromagnetic order along the backbone of a percolating cluster at p c . Thus, our 
estimates for d/z should be non-zero at p c for large transverse field strength T > Tm, since we expect 
the dynamical exponent z to diverge only at pfm(T), where pfm(T) is the critical concentration at 
transverse field strength T. The observation of a nonzero slope at p c implies that p c is not a special 
concentration for T > Tm- 

Now let us investigate the criticality of z ( |l9|) for T < Tm- The value of Dv p obtained from the 
percolation theory is known as Du p ~ 2.57. As is shown in Fig. 12, the criticality of d/z seems to 



be indeed compatible with the relation (19). For T > Tm the phase boundary is no longer parallel 
to the T-axis and the quantum fluctuations are strong enough to destroy the ferromagnetic order 
on the backbone of a percolating cluster at p c . We also analyzed the criticality of z for T > Tm- 
We do not know, however, the critical value of pfm (T). Thus we tried to fit the data in the 
form \p — pfm (T) \~ a - Minimizing the deviaton from the fitting form, we obtain pfm (1-5) — 0.67 
and a ~ 2.36, and pfm (2-0) ~ 0.77 and a ~ 2.24. The values of a are not far from for T < Tm- 
However, in order to conform the same universality class for the two regions, we need a more precise 
investigation which will be reported elsewhere. 

§4. Discussion and Summary 

The nature of the quantum Griffiths phase of the diluted transverse field Ising ferromagnet in two 
dimensions has been studied. Due to the Griffiths-McCoy singularity, the local susceptibility xiocal 
diverges, if d/z is smaller than 1. Thus, as in the case of a quantum spin glaflsP'El the Griffiths- 
McCoy singularity causes a divergence of various static quantities such as the zero-frequency local 
susceptibilities. In the classical case the Griffiths singularity is not strong enough to produce such 
a dramatic effect. 

So far we assume that the clusters are well defined geometrically and we have used the probability 
function P (N) given by percolation theory. However, even in the geometrically connected cluster 
the correlation function is reduced due to the quantum fluctuations, and when T becomes large 
we have to work with the probability distribution of "physical clusters" which are smaller than the 
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geometrical ones. Although we do not know P {N p h ys ) as a function of T at this moment, generally 
we expect to find the following senario: 

For T < Fm, P (Nphys) — P(N), and the mean cluster size diverges at p = p c . On the other 
hand for T > Fm, the mean size of the physical cluster diverges at pfm (F) > p c - In Fig. 11, we 
find that d/z vanishes at p c for the values F =0.7 and 1.0, which are both smaller than r^. On 
the other hand, for F > Fm they vanish at a larger value of p, which is considered to be pfm(T)- 
For T > Fm it is difficult to estimate the exponent correctly. Because the phase boundary at the 
zero temperature is not yet known correctly. The question of whether the critical exponent of this 
physical cluster may be different from the value of the exponent in the regions F < Fm is still open. 

For r ~ r c quantum fluctuation are very strong. In this case the present analysis which makes 
use of the broad distribution of xiocal does not apply for the small systems which we are able to 
investigate. The correlation function along the Trotter axis is not well developed, and the relation 
£ T ~ exp (CN) is no longer valid. The nature of this region which is dominated by quantum fluc- 
tuations will be studied elsewhere. 
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Fig. 1. Phase diagram of the transverse field Ising model on the square lattice for the finite temperatures. The 
symbols (o) denote the results of the series expansions by Elliott and Wood, and the symbols (+) denote the values 
obtained from the conventional method. In the diluted system. The critical line shrinks as p — > p c . 
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Fig. 2. The Monte Carlo sweep dependence of the order parameter (M 2 ^ at /3 = 0.2, T = 0.1 starting from the 
complete ferromagnetic state. 
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Fig. 3. Estimation of the critical r at (3 = 10. The cross (To ^ 3.05) at 1/n 2 = is the critical point obtained from 
the continuous time cluster algorithm. This value is compatible with the value r c ~ 3.03 obtained by extrapolation. 
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Fig. 4. Scaled uniform susceptibility. In this case, T c = 3.07, v = 0.63, and rj = 0.04. These estimates agree well 
with the series expansion results. 
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Fig. 5. Schematic phase diagram of the site diluted system. 
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Fig. 6. Schematic picture of a local cluster which forms a rod along the the Trotter axis. 



15 



Log F 
1 

0.5 h 



beta=40 o 

beta=70 * 

beta=140 □ 

slope=-0.416 -- 




local 



Fig. 7. Log-Log plot of the integrated distribution of the single site susceptibility Xiocai at T = 0.7, p = 0.30. In 
this case, we evaluate the distribution using 100 different 20 x 20 configurations. The dotted line is a fit to the 
straight line region of the data and has the slope —0.416. 
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Fig. 8. Similar to Fig. 6 but for F — 1.0, p = 0.40. In this case, we evaluate the distribution using 100 different 
20 x 20 configurations. The slope of the straight line is —0.309. 
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Fig. 9. Similar to Fig. 6 but for r = 1.5, p = 0.50. In this case, we evaluate the distribution using 100 different 
20 x 20 configurations. The slope of the straight line is —0.512. 
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Fig. 10. Similar to Fig. 6 but for T = 2.0, p = 0.66. In this case, we evaluate the distribution using 100 different 
20 x 20 configurations. The slope of the straight line is -0.391. 
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Fig. 11. Plot of d/z obtained by fitting the integrated distribution F (xiocal) f° r various T. For F = 0.7 and 1.0, 
d/z vanish at p ~ p c , but, for T = 1.5 and 2.0, these vanish at p > p c . When z > d, the average susceptibility will 
diverge. 
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Fig. 12. Extraction of the exponent from the data for d/z for T — 0.7 and 1.0, in which, Ap implies \p — p c \ 
is 0.59, the known percolation threshold). The value of Dv p obtained from the percolation theory is known 
Dv p ~ 2.57, and the criticality of d/z seems to be indeed compatible with the relation (21). 
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